clc; clear; close all;

% 参数设置
N_iter = 1000;
N_discard = 800;
m = 1;
c = 11.5;

num_a = 40;
num_b2 = 40;

a_list = linspace(0, 5, num_a);
b2_list = linspace(0, 2*pi, num_b2);

SE_mat = zeros(num_a, num_b2);

% 固定初始值
x0 = 0.3;
y0 = 0.5;
z0 = 0.6;

parfor ia = 1:num_a
    a = a_list(ia);
    SE_tmp = zeros(1, num_b2);
    for ib = 1:num_b2
        b2 = b2_list(ib);
        x = x0; y = y0; z = z0;

        X = zeros(N_iter,1);

        for k = 1:N_iter
            state = SIMM(x, y, z, a, b2, c);
            x = state(1);
            y = state(2);
            z = state(3);
            X(k) = x;
        end

        data = X(N_discard+1:end);
        r = 0.2 * std(data);
        SE_tmp(ib) = SE(data, m, r);
    end
    SE_mat(ia, :) = SE_tmp;
end

SE_min = min(SE_mat(:));
SE_max = max(SE_mat(:));
SE_norm = (SE_mat - SE_min) / (SE_max - SE_min);

% 构造网格坐标
[B2_grid, A_grid] = meshgrid(b2_list, a_list);

% 绘制三维曲面图
figure;
surf(B_grid, R_grid, SE_norm', 'EdgeColor', 'none');
set(gca, 'YDir', 'reverse'); 
xlabel('b_2');
ylabel('a');
zlabel('SE');
title('3D-SIMM系统的SE复杂度');
colorbar;
grid on;

